Stan code for the analysis of my Social and Cultural Dynamics exam

Plate notation

WAIC calculations

add_predictive_posterior <- function(dataset, model) {
    post <- extract.samples(model)
    dataset %>%
        mutate(theta_ind = purrr::pmap(list(n_ind, group, intensity, vowel), function(n,gro,int,vow) {
            iA = post$ia + post$ia_vowel[,vow] + post$ia_group[,gro]
            return(inv_logit(iA))
        }), theta_gro = purrr::pmap(list(n_gro, group, intensity, vowel), function(n,gro,int,vow) {
            gA = post$ga + post$ga_vowel[,vow] + post$ga_group[,gro]
            return(inv_logit(gA))
        }))
}

my_waic <- function(pred_post) {
    pred_post %>%
        mutate(
            lik_gro = purrr::pmap(list(k_gro, n_gro, theta_gro), function(k,n,theta){
                dbinom(k,n,theta)}), 
            lik_ind = purrr::pmap(list(k_ind, n_ind, theta_ind), function(k,n,theta){
                dbinom(k,n,theta)}), 
            lppd = purrr::map2_dbl(lik_gro, lik_ind, function(lik1, lik2) {
                log(mean(c(lik1, lik2)))}), 
            pwaic = purrr::map2_dbl(lik_gro, lik_ind, ~var2(c(.x, .y)))) %>%
        summarise(
            WAIC = -2 * (sum(lppd) - sum(pwaic)),
            WAIC_SE = sqrt(n() * 2 * var2(-2 * (lppd - pwaic))),
            pWAIC = sum(pwaic))
}


dataset = data_frame(
    n_gro = full.data$n_gro, 
    n_ind = full.data$n_ind, 
    k_gro = full.data$k_gro,
    k_ind = full.data$k_ind,
    group = full.data$group, 
    intensity = full.data$intensity, 
    vowel = full.data$vowel)

do_waic <- function(model, data = dataset) {
    add_predictive_posterior(data, model) %>%
        my_waic()
}
load("models/intercept.model")
load("models/conf_align_local.model")
load("models/indiscriminate_align_local.model")
load("models/conf_align_global.model")
load("models/cosine.model")
load("models/self_l.model")
load("models/align_l.model")
load("models/synergy_l.model")
load("models/self_entr.model")
load("models/align_entr.model")
load("models/synergy_entr.model")
load("models/self_l_and_local_confidence.model")
load("models/self_l_and_cosine.model")
load("models/self_l_cosine_local_confidence.model")
load("models/cosine_local_confidence.model")
load("models/self_l_x_local_confidence.model")
load("models/self_l_x_cosine.model")
load("models/cosine_x_local_confidence.model")

load("models/self_l_x_cosine_conf_local_full.model")
load("models/self_l_x_cosine_conf_local.model")
load("models/cosine_x_local_confidence_2.model")
load("models/self_l_x_cosine2.model")
load("models/self_l_x_local_confidence2.model")
load("models/self_l_x_local_confidence_only.model")


waic_list <- tribble(
    ~model,                             ~waic,
    "intercept",                        do_waic(intercept),
    "local confidence alignment",       do_waic(conf_align_local),
    "local indiscriminate alignment",   do_waic(indiscriminate_align_local),
    "global confidence alignment",      do_waic(conf_align_global),
    "word set cosine similarity",       do_waic(cosine),
    "self L",                           do_waic(self_l),
    "alignment L",                      do_waic(align_l),
    "synergy L",                        do_waic(synergy_l),
    "self entropy",                     do_waic(self_entr),
    "alignment entropy",                do_waic(align_entr),
    "synergy entropy",                  do_waic(synergy_entr),
    "self L + local confidence",        do_waic(self_l_and_local_confidence),
    "self L + cosine",                  do_waic(self_l_and_cosine),
    "self L + cosine + local conf",     do_waic(self_l_cosine_local_confidence),
    "cosine + local confidence",        do_waic(cosine_local_confidence),
    "self L * conf + cosine",           do_waic(self_l_x_local_confidence),
    "self L * cosine + local conf",     do_waic(self_l_x_cosine),
    "cosine * local confidence + self L", do_waic(cosine_x_local_confidence),
    "self L * local + self L * cosine", do_waic(self_l_x_local_confidence2),
    "self L * cosine + local * cosine", do_waic(self_l_x_cosine2),
    "self L * local + cosine * local",  do_waic(cosine_x_local_confidence_2),
    "self L * cosine * local 2way only",do_waic(self_l_x_cosine_conf_local),
    "self L * cosine * local",          do_waic(self_l_x_cosine_conf_local_full),
    "self L * local",                   do_waic(self_l_x_local_confidence_only)
    
) %>% unnest() %>%
    arrange(WAIC) %>%
    mutate(dWAIC = WAIC - WAIC[1],
           weight = exp(-0.5 * dWAIC),
           weight = weight / sum(weight)
    ) %>%
    # round the table to 3 digits
    mutate_at(vars(-model), function(co) {round(co, 3)})
waic_list
## # A tibble: 24 × 6
##                                 model     WAIC WAIC_SE pWAIC  dWAIC weight
##                                 <chr>    <dbl>   <dbl> <dbl>  <dbl>  <dbl>
## 1    self L * local + self L * cosine 1484.094 112.356 4.906  0.000  0.394
## 2              self L * conf + cosine 1485.220 112.065 4.859  1.126  0.224
## 3  cosine * local confidence + self L 1486.051 111.968 4.842  1.957  0.148
## 4             self L * cosine * local 1486.475 112.913 4.991  2.381  0.120
## 5    self L * cosine + local * cosine 1487.380 112.981 4.890  3.286  0.076
## 6        self L * cosine + local conf 1490.177 112.812 4.847  6.084  0.019
## 7     self L * local + cosine * local 1490.733 112.996 4.894  6.639  0.014
## 8        self L + cosine + local conf 1494.992 113.165 4.766 10.898  0.002
## 9   self L * cosine * local 2way only 1495.961 114.305 4.937 11.867  0.001
## 10                     self L * local 1496.377 113.061 4.792 12.283  0.001
## # ... with 14 more rows
# save("waic_list", file="models/waic.Rdata")

pred.post.data <- extract.samples(self_l_x_cosine_conf_local_full) %>%
    as.data.frame() %>%
    as_data_frame() %>%
    select(col_a : col_b7)

pred.post.data %>%
    gather() %>%
    # group_by()
    ggplot(aes(value, fill=key)) +
    geom_density(alpha=.4, colour=NA) +
    facet_wrap(~ key, scales="free")

pairs(pred.post.data)

# confidence alignment back-scaling
# m = mean(text$local_alignment_confidence)
# s = sd(text$local_alignment_confidence)
# 
# data_frame(local_confidence = seq_range(full.data$local_confidence, 10)) %>%
#     mutate(samples = rep(list(pred.post.data), n()),
#            pred.post = purrr::map2(local_confidence, samples, ~.y$col_a + .y$col_b_confidence * .x),
#            local_confidence = (local_confidence * s) + m) %>%
#     unnest(pred.post) %>%
#     ggplot(aes(local_confidence, pred.post)) +
#     geom_line(stat="summary", fun.y=mean) +
#     geom_ribbon(stat="summary", fun.data=mean_hpdi, alpha=.3) +
#     geom_point(aes(local_confidence, y = .97), data = data.frame(local_confidence = (full.data$local_confidence * s) + m),
#                shape="|", size=3) +
#     theme_tufte() +
#     theme(axis.ticks = element_blank()) +
#     labs(x = "Local confidence alignment: Proportion of confidence expressions aligned with previous sentence",
#          y = "Predicted collective benefit, mean +- 89% HPDI",
#          title = "Counterfactual predictive posterior for explicit metacognition (confidence alignment)",
#          subtitle = "Actual dyads' confidence alignment as x-axis") +
#     scale_y_continuous(breaks = seq(.96,1.1, .02))

# ggsave("plots/pred_post_local_confidence.png")
posterior = extract.samples(self_l_x_cosine_conf_local_full)

pred <- expand.grid(local_alignment_confidence = seq(-2,2,length.out = 10),
          self.l = seq(-2,2, 2),
          cosine_word_set = seq(-2,2,2)) %>%
    mutate(collective_benefit = purrr::pmap(list(local_alignment_confidence, self.l, cosine_word_set), function(loc, sel, cosi) {
        posterior$col_a +
            posterior$col_b1 * sel +
            posterior$col_b2 * loc +
            posterior$col_b3 * cosi + 
            posterior$col_b4 * loc * sel +
            posterior$col_b5 * cosi * sel +
            posterior$col_b6 * cosi * loc +
            posterior$col_b7 * cosi * loc * sel
    })) %>%
    plyr::adply(1, function(dat) {
        # print(dat)
        mean_hpdi(dat$collective_benefit[[1]])
    }) %>%
    select(-collective_benefit)



pred %>%
    ggplot(aes(local_alignment_confidence, ymin=ymin, y=y, ymax=ymax)) +
    geom_line() +
    geom_ribbon(alpha=.3) +
    theme_tufte() +
    theme(axis.ticks = element_blank()) +
    labs(x = "Local confidence alignment: Proportion of confidence expressions aligned with previous sentence (z-scale)",
         y = "Predicted collective benefit, mean +- 89% HPDI",
         title = "Counterfactual predictive posterior for the 3 way interaction model",
         subtitle = "Facets are z-scaled self L (chat predicability) and cos distance between word sets of interaction partners") +
    # scale_y_continuous(limits = c(0,1)) +
    coord_cartesian(ylim = c(0,6)) +
    facet_grid(cosine_word_set ~ self.l, labeller=label_both)

# ggsave("plots/3way_interaction.png")
posterior = extract.samples(self_l_x_local_confidence2)

pred <- expand.grid(local_alignment_confidence = seq(-2,2,length.out = 10),
          self.l = seq(-2,2, 2),
          cosine_word_set = seq(-2,2,2)) %>%
    mutate(collective_benefit = purrr::pmap(list(local_alignment_confidence, self.l, cosine_word_set), function(loc, sel, cosi) {
        posterior$col_a +
            posterior$col_b1 * sel +
            posterior$col_b2 * loc +
            posterior$col_b3 * cosi + 
            posterior$col_b4 * loc * sel +
            posterior$col_b5 * cosi * sel;
    })) %>%
    plyr::adply(1, function(dat) {
        # print(dat)
        mean_hpdi(dat$collective_benefit[[1]])
    }) %>%
    select(-collective_benefit)



pred %>%
    ggplot(aes(local_alignment_confidence, ymin=ymin, y=y, ymax=ymax)) +
    geom_line() +
    geom_ribbon(alpha=.3) +
    theme_tufte() +
    theme(axis.ticks = element_blank()) +
    labs(x = "Local confidence alignment: Proportion of confidence expressions aligned with previous sentence (z-scale)",
         y = "Predicted collective benefit, mean +- 89% HPDI",
         title = "Counterfactual predictive posterior for: self L * local + self L * cosine",
         subtitle = "Facets are z-scaled self L (chat predicability) and cos distance between word sets of interaction partners") +
    # scale_y_continuous(limits = c(0,1)) +
    coord_cartesian(ylim = c(0,3)) +
    facet_grid(self.l ~ cosine_word_set, labeller=label_both)

posterior = extract.samples(self_l_x_local_confidence)

pred <- expand.grid(
        local_alignment_confidence = seq(-2,2,length.out = 10),
        self.l = seq(-2,2, 2),
        # cosine_word_set = seq(-2,2,length.out=10)
        # self.l = 0,
        # local_alignment_confidence = 0,
        cosine_word_set = 0
        ) %>%
    mutate(collective_benefit = purrr::pmap(list(local_alignment_confidence, self.l, cosine_word_set), function(loc, sel, cosi) {
        posterior$col_a +
            posterior$col_b1 * sel +
            posterior$col_b2 * loc +
            posterior$col_b3 * cosi + 
            posterior$col_b4 * loc * sel;
    })) %>%
    plyr::adply(1, function(dat) {
        # print(dat)
        mean_hpdi(dat$collective_benefit[[1]])
    }) %>%
    select(-collective_benefit)


pred %>%
    ggplot(aes(local_alignment_confidence, ymin=ymin, y=y, ymax=ymax)) +
    geom_line() +
    geom_ribbon(alpha=.3) +
    theme_tufte() +
    theme(axis.ticks = element_blank()) +
    labs(x = "Local confidence alignment: Proportion of confidence expressions aligned with previous sentence (z-scale)",
         y = "Predicted collective benefit, mean +- 89% HPDI",
         title = "Counterfactual predictive posterior for: self L * local + cosine",
         subtitle = "Facets are z-scaled self L (chat predicability) and cos distance between word sets of interaction partners") +
    # scale_y_continuous(limits = c(0,1)) +
    coord_cartesian(ylim = c(0.7,1.4)) +
    facet_wrap(~ self.l, labeller=label_both)

# gA = ga + ga_vowel[vowel] + ga_group[group]; 
# gB = gb + gb_vowel[vowel] + gb_group[group];
# gtheta[i] =  gA[i] + gB[i] * intensity[i];

# expand.grid(intensity = unique(full.data$intensity),
#             vowel = 1:4,
#             group = 1:9) %>%
#     mutate(Group = pmap(list(intensity, group, vowel), function(i,g,v){
#         post$ga + post$ga_vowel[,v] + post$ga_group[,g] +
#             (post$gb + post$gb_vowel[,v] + post$gb_group[,g]) * i
#     }), Individual = pmap(list(intensity, group, vowel), function(i,g,v){
#         post$ia + post$ia_vowel[,v] + post$ia_group[,g] +
#             (post$ib + post$ib_vowel[,v] + post$ib_group[,g]) * i
#     })) %>%
#     left_join(data.frame(group = 1:9, confidence = full.data$local_confidence)) %>%
#     gather("who", "theta", c(Group, Individual)) %>%
#     mutate(n = 10,
#            theta = purrr::map(theta, ~inv_logit(.))) %>%
#     mutate(vowel = c("e", "i", "o", "u")[vowel]) %>%
#     unnest(theta) %>%
#     ggplot(aes(intensity, theta, group=str_c(group, who), colour=confidence)) +
#     # geom_ribbon(aes(group=NA), stat="summary", fun.data=mean_hpdi, alpha=.3) +
#     geom_line(stat="summary", fun.y=mean, alpha=.7) +
#     geom_line(aes(group=NA), stat="summary", fun.y=mean, colour="red") +
#     facet_wrap(~ vowel, labeller = label_both) +
#     theme_few() +
#     labs(x = "Stimulus intensity (dB)", y = "Proportion correctly identified")

comparing the four models from optimally interacting ..

Coin flip model

(s1 + s2) / 2

Behaviour and feedback model

max(s1, s2)

weighted confidence sharing

(s1 + s2) / (2 ^ 1/2)

direct signal sharing

(s1 ^ 2 + s2 ^ 2) ^ 1/2